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Accessing phonon polaritons in hyperbolic crystals by ARPES 
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Recently studied hyperbolic materials host unique phonon-polariton (PP) modes. The ultra- 
short wavelengths of these modes, which can be much smaller than those of conventional exciton- 
polaritons, are of high interest for extreme sub-diffraction nanophotonics schemes. Polar hyperbolic 
materials such as hexagonal boron nitride can be used to realize strong long-range coupling between 
PP modes and extraneous charge degrees of freedom. The latter, in turn, can be used to control 
and probe PP modes. Of special interest is coupling between PP modes and plasmons in an adja¬ 
cent graphene sheet, which opens the door to accessing PP modes by angle-resolved photoemission 
spectroscopy (ARPES). A rich structure in the graphene ARPES spectrum due to PP modes is 
predicted, providing a new probe of PP modes and their coupling to graphene plasmons. 


Introduction .—The intrinsic hyperbolic character [1] of 
hexagonal boron nitride (hBN) grants a unique platform 
for realizing deep-subwavelength nanophotonic schemes. 
Key to these developments are phonon-polariton (PP) 
modes that exist within reststrahlen frequency bands [2, 
3], characterized by wavelengths that can be as small as 
1-100 nm. Highly directional, these modes exhibit deep 
sub-diffraction confinement of light with wavelengths far 
shorter than those of exciton-polaritons in semiconduc¬ 
tor microcavities [4]. PPs have been shown to propagate 
with low losses [2, 3] besting artificial metallic-resonator 
metamaterial schemes, and opening the door to hyper- 
lensing [5, 6]. 

Key to harnessing PP modes is gaining access to their 
response over a wide wavenumber and energy band- 
widths. However, to date PPs have only been studied 
within a small frequency range limited by laser choice 
(e.g. 167 meV ^ huj < 198 meV via scattering- 

type near-field optical spectroscopy technique [2]), or 
at specific wavelengths fixed by the sample geometry 
via Fourier transform infrared spectroscopy of nanofabri- 
cated nanopillars [3]. New approaches allowing to resolve 
the PP modes at shorter wavelengths and over a broad 
range of energies are therefore highly desirable. 

Here we describe an angle-resolved photoemission 
spectroscopy (ARPES) [9] scheme to achieve broadband 
energy-resolved access to ultra-short wavelength PPs in 
hBN. At first glance, ARPES access to PPs in a wide- 
bandgap insulator (hBN) where no free carriers are avail¬ 
able may seem counterintuitive. However, the key to 
our protocol lies in coupling PPs to charge degrees of 
freedom in a conductor (e.g. graphene) placed nearby 
the hyperbolic crystal of interest (hBN), prepared in a 
slab geometry. Strong coupling [10-15] between hBN 
Fabry-Perot PP modes and the collective charge oscilla¬ 
tions (i.e. Dirac plasmons [16]) in a doped graphene sheet 
placed over a hBN slab gives rise to new channels for 
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FIG. 1. (Color online) Signatures of PP modes in the quasi¬ 
particle spectral function A{k, cj) of a doped graphene sheet 
placed over a hBN slab, obtained from Eqs. (3) and (4). Note 
the black linearly-dispersing quasiparticle bands, which dis¬ 
play a clear Dirac crossing labeled by (1), and the broad spec¬ 
tral feature labeled by (2) due to the emission of the plasmon- 
phonon polariton mode with highest energy in Fig. 2(a). The 
Fermi energy is positioned at Ct; = 0. Emission of polariton 
modes [see Fig. 2(a)] by the holes created by photo-excited 
electrons gives rise to four dispersive satellite bands running 
parallel to the main quasiparticle bands (marked by red ar¬ 
rows). The feature (2) is mainly plasmonic,whereas the satel¬ 
lite bands, crossing at /c = 0 between features (1) and (2), are 
entirely due to Fabry-Perot hBN phonon-polariton modes. 
Parameters used: Fermi energy ef = 400 meV, hBN slab 
thickness d = 60 nm, Ca = 1 (vacuum), eb = 3.9 (Si02). The 
colorbar refers to the values of hA{k,uj) in eV. 
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quasiparticle decay yielding a rich structure of dispersive 
satellite features—marked by red arrows in Fig. 1— in the 
graphene ARPES spectrum A{k,uj). Since hBN Fabry- 
Perot PPs are controlled by slab thickness, the composite 
G/hBN structure features a novel ARPES spectrum with 
features that are highly tunable by thickness. 

The greatest practical advantage of this approach is 
that ARPES achieves extreme resolution over a wide 
range of wave vectors k (from the corners K, K' of the 
graphene Brillouin zone to the Fermi wave number 
in graphene) and energies huj^ with all energies below 
the Fermi energy being probed simultaneously. This 
gives an additonal benefit, besides tunability, in that 
the entire range of frequencies and wavenumbers can 
be covered within a single experiment. It is remarkable 
that a one-atom-thick conducting material like graphene, 
once placed over an insulating hyperbolic crystal, enables 
ARPES studies of PP modes over the full range of wave 
vectors and energies of interest. 

From a more fundamental perspective, ARPES will 
also be an ideal tool to investigate whether effective 
electron-electron interactions mediated by the exchange 
of PPs are capable of driving electronic systems towards 
correlated states. Finally, looking at our results from the 
point of view of graphene optoelectronics, one can en¬ 
vision situations in which the tunable coupling between 
graphene quasiparticles and the complex excitations of 
its supporting substrate can be used to achieve control 
over the spectral properties of graphene carriers, includ¬ 
ing their decay rates, renormalized velocities, etc. This 
degree of tunability may have important implications on 
the performance of graphene-based photodetectors [17]. 
Phonon and plasmon-phonon polaritons. —We consider 
a vertical heterostructure—see inset in Fig. 2(b)— 
composed of a graphene sheet located at 2 ; = 0 and placed 
over a homogeneous anisotropic insulator of thickness d 
with dielectric tensor e = diag(ea,, e^, 6;^). Homogeneous 
and isotropic insulators with dielectric constants Ca and 
Cb fill the two half-spaces 2 ; > 0 and z < —d, respec¬ 
tively. The Fourier transform Vq^uj of tho Coulomb inter¬ 
action potential, as dressed by the presence of a uniaxial 
{^y = ^x) dielectric, is given by 

y _ _ + €],ta.nh{qdy'ex/€z) _ 

+ ebCa) tanh{qd^yex/€z)/{2€) ’ 

( 1 ) 

where Vq = 27re^/(ge) with e = (ea + eb)/2 is the ordi- 
nary 2D Coulomb interaction potential. A more general 
equation, which is also valid in the case ^ can be 
found in Sect. I of Ref. 18. 

In the case of hBN, the components and of 
the dielectric tensor have an important dependence on 
frequency uj in the mid infrared [19]. The simplest 
parametrization formulas for Cx^z = are reported 

in Sect. I of Ref. 18 and have been used for the numer¬ 
ical calculations. More realistic parametrizations can be 


found in the Supplementary Information of Ref. 13. 

Standing PP modes [2] correspond to poles of the 
dressed interaction Vq^uj inside the reststrahlen bands. 
These can be found by looking at the zeroes of the denom¬ 
inator in Eq. (1), ^J\ex{u^) ez{uJ)\ + {2e)~^[ex{uj)ez{u:) + 
CbCa] i^in[qd^J\ex{oJ)/ez{oj)\] = 0. Illustrative numerical 
results for d = 10 nm and d = 60 nm are reported in 
Fig. 1 of Ref. 18. Analytical expressions, which are valid 
for qd 1 and qd ^ are available [18] in the case 
in which phonon losses in hBN are neglected. For suffi¬ 
ciently thick hBN slabs, there can be modes with group 
velocity equal to the graphene Fermi velocity up- 

Standing PP modes in a hBN slab couple to Dirac plas- 
mons in a nearby graphene sheet. Such coupling is cap¬ 
tured by the random phase approximation (RPA) [20]. 
In the RPA, one introduces the dynamically screened in¬ 
teraction 

w = _ (2) 

e{q,u)- 1-Vq,^xo{q,coy 


Here e{q^uj) is the RPA dielectric function and Xo{q^^) 
is the density-density response function of a 2D massless 
Dirac fermion fluid [21]. While the poles of Vq^uj physi¬ 
cally yield slab PP modes, new poles of Wq^uj emerge from 
electron-phonon interactions. These are weakly-damped 
solutions uj = Clq — i0+ of the equation e{q^u) = 0. 
We have solved this equation numerically and illustra¬ 
tive results for = 400 meV and d = 60 nm are shown 
in Fig. 2(a). (Results for different values of and d 
can be found in Sect. HI of Ref. 18.) Solid lines repre¬ 
sent plasmon-phonon polaritons that emerge from the hy¬ 
bridization between the Dirac plasmon [16] in graphene 
(dashed line) and standing PP waves in the hBN slab. 
The solid red lines denote three polariton branches with 
a strong degree of plasmon-phonon hybridization. On 
the contrary, black solid lines denote practically unhy¬ 
bridized slab PP modes. We clearly see that there are 
several plasmon-phonon polariton modes (green circles) 
with group velocity equal to up. These modes couple 
strongly to quasiparticles in graphene, as we now pro¬ 
ceed to demonstrate. 

Quasiparticle decay rates. —An excited quasiparticle with 
momentum k and energy huj^ created in graphene in 
an ARPES experiment [22-26], can decay by scattering 
against the excitations of the Fermi sea, i.e. electron-hole 
pairs and collective modes. The decay rate h/Tx{k, uj) for 
these processes can be calculated [20] from the imaginary 
part of the retarded quasiparticle self-energy Tix{k^uj), 
i.e. h/Tx{k,uj) = —2lm[T^x{k,uj)]. In the RPA and at 
zero temperature we have [27, 28] 


Im [I]a(A:,w)] = 


\f ^ 


d?q 


Im 






AA' 


X [Q{hw - ix',k+q) - 0(-CA',fc+q)] • (3) 


Here Pxx' = [1 + AA'cos (^fc^fc+q)]/2 is the chirality fac¬ 
tor [27, 28], = Afiup/c — sf is the Dirac band energy 
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measured from the Fermi energy (A, A' = ±1), and 
&{x) is the usual Heaviside step function. The quantity 
fuj is also measured from the Fermi energy and, finally, 
^k,k-\-q is the angle between k and k-\-q. Eq. (3) reduces 
to the standard Fermi golden rule when only terms of 
0{Vq^^) are retained. Physically, it describes the decay 
rate of a process in which an initial state with momen¬ 
tum k and energy hw (measured from 5 f) decays into 
a final state with momentum k q and energy 
(measured from 5 f). For cj < 0, the self-energy expresses 
the decay of holes created inside the Fermi sea, which 
scatter to a final state, by exciting the Fermi sea. Fermi 
statistics requires the final state to be occupied so both 
band indices A' = ±1 are allowed in the case 5 f > 0 that 
we consider here. Since ARPES measures the properties 
of holes produced in the Eermi sea by photo-ejection, 
only cj < 0 is relevant for this experimental probe in an 
n-doped graphene sheet. 

It is convenient to discuss the main physical features 
of Im [Ea(/c, cj)] for an initial hole state with momen¬ 
tum Aj = 0. In this case, the 2D integral in Eq. (3) 
reduces to a simple ID quadrature. The initial hole 
energy is = hw Sy- The final hole energy is 
Ei = ^x',q + ^F = \'hvYq> When the difference = 

Ei — Ei is equal to the real part of the mode energy 
hQq^ the initial hole, which has been left behind after 
the photo-ejection of an electron, can decay by emit¬ 
ting a plasmon-phonon polariton. Since hClq > hv^q, 
but Ax'^q < hv^q for intraband transitions, an initial 
hole state with Ei < 0 (i.e. initial hole state in valence 
band) can decay only into a final hole state with Ef > 0 
(i.e. final hole state in conduction band). In particular, 
when dQqjdq = h~^dAx',q/dq = ATf, such decay pro¬ 
cess is resonant When these conditions are met, the 
inter-band contribution to Im [I1 a(0,cc;)] peaks at a char¬ 
acteristic value of uj and the Kramers-Kronig transform 
Re[EA(0, cj)] changes sign rapidly around that frequency. 
Within RPA, a satellite quasiparticle emerges [29], which 
is composed by a hole that moves with the same speed 
of a plasmon-phonon polariton. This is a solution of the 
Dyson equation, distinct from the ordinary quasiparticle 
solution that becomes the Landau pole of the one-body 
Green’s function as k ^ and u; ^ 0. 

The quantity Im [Ea(0,ci;)], calculated from Eq. (3), is 
plotted as a function of uj in Eig. 2(b), for sf = 400 meV 
and d = 60 nm. (The dependence of the decay rate 
on sy and d is discussed in Sect. Ill of Ref. 18.) We 
clearly see several peaks in Im [Ea( 0, cj)] for huu < —sy 
(Ei < 0), which occur at values of huj that are in a 
one-to-one correspondence with the “resonant” plasmon- 
phonon polaritons, i.e. polaritons with group velocity 
equal to vy, shown in Eig. 2(a). Indeed, as stated above, 
peaks in Im [Ea( 0, cj)] are expected at values of huj — 
marked by green vertical lines in Eig. 2(b) —given by 
hjjj = hvYq^ — Sy — hQq*, where q'^ is the wave num¬ 
ber at which the resonance condition dtlqjdq = vy is 
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FIG. 2. (Color online) Panel (a) Dispersion relation Qq of 
hybrid plasmon-phonon polaritons (solid lines) with param¬ 
eters as in Fig. 1. The dashed line represents the dispersion 
relation of a Dirac plasmon [16] in graphene, in the absence 
of hBN phonons. Horizontal cyan areas denote the hBN rest- 
strahlen bands. The grey-shaded area represents the intra¬ 
band particle-hole continuum in graphene. Green filled cir¬ 
cles represent the points where the plasmon-phonon polariton 
group velocity equals the graphene Fermi velocity i’f- Panel 
(b) The quantity —Im[EA(fc, cj)] (in units of ef and evaluated 
at = 0 ) is shown as a function of the rescaled frequency 
hw/eF- Green vertical lines denote the values of hujeF at 
which a plasmon-phonon polariton peak is expected. The ver¬ 
tical axis is in logarithmic scale. The inset shows a side view of 
the vertical heterostructure analyzed in this work. Panel (c) 
Same as in panel (b) but in the absence of dynamical screening 
due to electron-electron interactions in graphene: these nu¬ 
merical results have been obtained by replacing Wq^uj Fq,u; 
in Eq. (3). A polaron peak is clearly visible. 


satisfied. Eor example, the resonant mode at highest en¬ 
ergy in Eig. 2(a), which occurs at q'^ ~ 0.26 nm“^ and 
energy ^ 0.36 eV, yields a peak in Im [Ea(0,co’)] at 
hjjjeY ~ —1.5, see Eig. 2(b). 

Gomparing Eig. 2(b) with Eig. 2(c), we clearly see the 
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FIG. 3. (Color online) Panel (a) The quasiparticle spectral 
function A(k^uj) evaluated at = 0 (|fc| = 10“^ /cf has been 
used in the numerical calculations) is plotted as a function of 
hw. This plot refers to d = 10 nm. The other parameters are 
as in Fig. 1. Panels (b) and (c) Dependence on the hBN slab 
thickness d of the spectral function features highlighted by 
vertical orange-shaded regions in panel (a). Different curves 
correspond to values of d on a uniform mesh from d = 10 nm 
to d = 60 nm. Arrows indicate how spectral features evolve 
by increasing d. 


role of dynamical screening due to electron-electron inter¬ 
actions in graphene. For e{k^uj) = 1, the off-shell decay 
rate Im [Sy(0,ci;)] shows only a polaron peak, due to the 
emission of a Fabry-Perot PP mode with group velocity 
equal to see Fig. 1 in Ref. 18. 

At k 7 ^ 0, the conduction and valence band 

Im[SA(A5, cj)] plasmon-phonon polariton peaks broaden 
and separate, because of [27] the dependence on scatter¬ 
ing angle of ^x',k-\-q and the chirality factor Tw', which 
emphasizes k and q in nearly parallel directions for con¬ 
duction band states and k and q in nearly opposite direc¬ 
tions for valence band states. As a result, the conduction 
band plasmon-phonon polariton peak moves up in energy 


while the valence band peak moves down. 

Quasiparticle spectral function .—An ARPES exper¬ 
iment [9] probes the quasiparticle spectral function 
A{k,u}) = -7r“i y;;^Im[G'A(fe,w)] = Ea=±i w) 

of the occupied states below the Fermi energy. Here 
Gx{k,uj) is the one-body Green’s function in the band 
representation and 

Ax =_ _ 

^ - fx,k/h - ReJ^x/h)^ + (ImJ^x/h)^ ' ^ ^ 

In writing Eq. (4) we have dropped explicit reference 
to the k^uj variables. The real part Re[EA(A5, u;)] of 
the quasiparticle self-energy can be calculated, at least 
in principle, from the Kramers-Kronig transform of 
Im[EA(/c, cj)]. A more convenient way to handle the 
numerical evaluation of Re[EA(/c, cj)] is to employ the 
Quinn-Eerrell line-residue decomposition [30]. 

Our main results for the quasiparticle spectral func¬ 
tion A{k,uj) of a doped graphene sheet placed on a hBN 
slab are summarized in Eig. 1 and Eig. 3. We clearly 
see that the presence of the hBN substrate is respon¬ 
sible for the appearance of a family of sharp dispersive 
satellite features associated with the presence of PPs and 
plasmon-phonon polaritons. This is particularly clear in 
the one-dimensional cut at /c = 0 of A{k,uj) displayed 
in Eig. 3(a) for d = 10 nm. All the sharp structures 
between the ordinary quasiparticle peak slightly below 
huj = —0.4 eV and the peak at huj ~ —0.7 eV, which 
is mostly plasmonic in nature, are sensitive to the de¬ 
tailed distribution and dispersion of Eabry-Perot PP in 
the hBN slab, and therefore to the slab thickness d. This 
is clearly shown in Eig. 3(b) and (c), where we see shifts 
of these peaks of several meV, when d is changed from 
d =10 nm to d = 60 nm, while keeping constant. 

In summary, we have studied the coupling between 
standing phonon-polariton modes in a hyperbolic crys¬ 
tal slab and the plasmons of the two-dimensional mass¬ 
less Dirac fermion liquid in a nearby graphene sheet. 
We have shown that this coupling yields a complex 
spectrum of (plasmon-phonon) polaritons, see Eig. 2(a). 
Plasmon-phonon polaritons with group velocity equal 
to the graphene Eermi velocity couple strongly with 
graphene quasiparticles, enabling ARPES access to PP 
modes in hyperbolic crystal slabs, as shown in Eigs. 1 
and 3. Recent progress [31] in the chemical vapor de¬ 
position growth of large-area graphene/hBN stacks on 
Cu(lll) in ultrahigh vacuum and the ARPES character¬ 
ization of the resulting samples makes us very confident 
on the observability of our predictions. Our findings sug¬ 
gest that appropriate coupling of graphene to substrates 
which allow strong plasmon-phonon hybridization could 
open the route to the manipulation of carriers’ spectral 
properties, paving the way for novel device functionali¬ 
ties. 
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In this Supplemental Material File we present i) analytical details on the Coulomb interaction 
potential dressed by the presence of an anisotropic dielectric, ii) numerical and analytical results on 
the dispersion of the phonon-polariton modes, and iii) numerical results on the quasiparticle decay 
rates and spectral function as the Fermi energy and the dielectric thickness are varied. 


I. ELECTROSTATICS AND 
FREQUENCY-DEPENDENCE OF THE HBN 
DIELECTRIC TENSOR 

We consider a vertical heterostructure [inset of 
Fig. 2(b) in the main text] composed of a graphene sheet 
located at 2 ; = 0 and placed over a homogeneous but 
anisotropic insulator of thickness d with dielectric tensor 
€ = dmg{ex^Cy,ez). Homogeneous and isotropic insula¬ 
tors with dielectric constants Ca and Cb fill the two half¬ 
spaces 2 ; > 0 and z < —d, respectively. We calculate the 
electrical potential created by an electron in graphene for 
arbitary values of Cx 7 ^ Cy and then specialize our result 
to the case of a uniaxial crystal with Cx = Cy, such as 
hBN. 

The electrical potential can be calculated from the 
following elementary approach. The displacement field 
D{r,z) must satisty the condition 'V 'D{r,z) =0 every¬ 
where in space. However, the presence of an electron with 
charge density —eS^{r)S{z) at z = 0 implies a discontinu¬ 
ity of the normal component of the displacement field 
across 2 ; = 0 , while the tangential components Ex^Ey of 
the electric field E{r,z) must be continuous. 

Since the electric field E{r,z) is irrotational every¬ 
where in space, we can introduce the electric potential 
U(r, z) in the three regions of space 2 : > 0 , —d < 2 : < 0 , 
and 2 : < —d. The Laplace equation —exd‘^V{r,z) — 
€ydlV{r,z)-€,dlV(r,z) = 0 in the anisotropic dielectric 
(i.e. for —d < 2 ; < 0 ) can be reduced [ 1 ] to an ordinary 
Laplace equation by scaling x xjs/ef, y yj 
and 2 ; ^ zj Imposing the aforementioned boundary 
conditions and carrying out elementary algebraic steps, 
we find the 2D Fourier transform Vq of the electrical po¬ 
tential: 

y ^ + ebg^ tanh (cIk) _ 

e^qK + (eaCbg^ + e^K^) tanh {dK)/{2e) 

Here, q = |q|, q = {qx,qy), k = {e^qlhz + 
and ifq = —27Tel{qe) with e = (ea + eb)/2. Eq. (1) 
is the most important result of this Section and repro¬ 
duces all known elementary results. For example, in the 
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TABLE 1. Microscopic parameters [4] entering Eq. (3). 


limit d ^ 0 Eq. (1) yields the potential of an electron 
in a graphene sheet embedded in two homogeneous and 
isotropic dielectrics, i.e. Vq ipq. Similarly, for any 
finite d and for an isotropic medium with Cx = Cy = 
Cz = €[, it is easy to check that Eq. (1) yields [2] Vq 
—dTrecosh {tJq) cosh {qd + '0h)/[Q^i sinh(gd + 77a TT/b)] with 
2r/a,b = In [(ei -h ea,b)/(ei - ea,b)]- The ratio on the right- 
hand side of Eq. (1) can be interpreted as a form fac¬ 
tor due to the presence of the anisotropic insulator slab, 
which dresses the “bare” 2D Coulomb potential ^pq. Eor 
graphene on hBN, the relevant limit of Eq. (1) is that of 
a uniaxial dielectric medium with Cy = In this case 
Eq. (1) reduces to 

y _ _ ^/exez + ebtanh(gtiye3;/e^) _ 

+ {ix^z + EbCa) ta.nh{qd^/ZJe^)/ (2e) 

( 2 ) 

We have modified the notation of the Coulomb po¬ 
tential in Eq. (2) to explicitly indicate its dependence 
on frequency. Indeed, in the case of hBN, the compo¬ 
nents of the dielectric tensor have an important depen¬ 
dence on frequency in the mid infrared, which is usually 
parametrized in the following form [3] 

/ X _ , Q,0 — Q,c» /o\ 

with i = X 01 z. Here and e^^oo are the static and 
high-frequency dielectric constants, respectively, while 
ujJ is the transverse optical phonon frequency in the di¬ 
rection £. The longitudinal optical phonon frequency 
cjj" satisfies the Lyddane-Sachs-Teller relation = 
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FIG. 1. (Color online) Dispersion of phonon-polaritons 
modes in a hBN slab. Poles of the dressed electrical potential 
Vq,uj in Eq. (2) are shown as functions of the wave vector q. 
Here and in the following figures we use Ca = 1 (vacuum) and 
Ch = 3.9 (Si02). Panel (a) is for d = 10 nm, while panel (b) 
is for d = 60 nm. Shaded areas denote upper and lower rest- 
strahlen bands. The green filled circle represents the point 
where the group velocity of the mode equals the graphene 
Fermi velocity vf- The 10 nm-thick slab does not support 
modes with group velocity equal to r'p. 

ujJ \/e^^o/Q,cx)- The parameters in Eq. (3) are listed 
in Table I and have been taken from recent measure¬ 
ments [4] on high-quality bulk hBN [5] . A simple inspec¬ 
tion of Eq. (3) in the limit 7 ^ ^ 0 shows that, as we 
stated earlier, hBN is a hyperbolic material: in the lower 
(upper) reststrahlen band, which is defined by the in¬ 
equality ujJ < uj < uj^ (cjJ < uj < cj^), the quantities 
take negative values. 


FIG. 2. (Color online) Panel (a) Plasmon-phonon polari- 
ton dispersion for a doped graphene sheet (sf = 400 meV) 
on a 10 nm-thick hBN slab. The dashed line represents the 
dispersion relation of a Dirac plasmon in graphene in the ab¬ 
sence of hBN phonons. Horizontal shaded areas and green 
filled circles have the same meaning as in Fig. 1. The density 
plot shows the imaginary part of the non-interacting polar¬ 
ization function in graphene and the corresponding colorbar 
is in units of the density of states at the Fermi energy. Panel 
(b) The quantity —Im[EA(fc, cj)] (in units of sf and evaluated 
at = 0) is shown as a function of the rescaled frequency 
hw/eF- Green vertical lines denote the values of bwleF at 
which a plasmon-phonon polariton peak is expected. 

Here, we report the results for the modes in the upper 
reststrahlen band. Similar expressions can be obtained 
for the lower reststrahlen band. For gd ^ 1, we find that 
one mode approaches the bottom of the reststrahlen band 
linearly 


II. ANALYTICAL RESULTS ON THE 
DISPERSION PHONON-POLARITON MODES 

In the limit 7 ^ ^ 0 it is possible to obtain analyt¬ 
ical expressions for the dispersion of phonon-polariton 
modes in a hBN slab (without a graphene sheet placed 
on top of it). We remind the reader that these modes are 
the solutions uu = uuq of the equation y^\ex{uj)ez{uj)\ + 
{2e)~'^[€x{ui)€z{u) + eb€a]ta.n[qd^y\€JJJ)J€JJJ)\] = 0 . 
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while the other modes decrease quadratically 



^x,0 ^x,oo 



( 4 ) 


( 5 ) 


In the limit qd ^ 1 (with n = 1, 2,...), instead we find 
that all modes approach the upper end of the reststrahlen 
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FIG. 3. (Color online) Same as Fig. 2 with parameters 
= 150 meV and d = 60 nm. 



ep [eV] 


FIG. 4. (Color online) A 2D color plot of the quantity 
—Im[EA(0,cj)] as a function of HujIsf and ep- Data in this 
plot refer to a hBN slab of thickness d = 60 nm. As usual, 
Ca = 1 and Cb = 3.9. The horizontal dashed line marks the 
expected position of a plasmon-phonon polariton peak in the 
absence of the hBN slab. The vertical dotted lines mark the 
extremes of the upper reststrahlen band. 


band asymptotically 


iUa 





^x,0 ^x,oo 
^x,oo^x,0 



(6) 


III. ADDITIONAL NUMERICAL RESULTS 

The following figures contain additional numerical re¬ 
sults with respect to those included in the main text. 

Figs. 2 and 3 show the plasmon-phonon polariton mode 
dispersion and the frequency dependence of the imagi¬ 
nary part of the quasiparticle self-energy. With respect 
to Figs. 2(a) and 2(b) in the main text, different values 
of hBN thickness and graphene Fermi energy are used in 
Fig. 2 and 3, respectively. In Fig. 2, we see that the dis¬ 
persion of the plasmon-phonon polaritons changes with 
the hBN slab thicknes, and so does the wave vector where 
the group velocity of the modes equals the graphene 
Fermi velocity vp. However, the connection between the 
modes slope and the peaks in the imaginary part of the 
quasiparticle self-energy persists. In Fig. 3, the lower 
chemical potential induces a less steep dispersion of the 
bare plasmon mode. This reflects in the fact that all the 
plasmon-phonon polariton modes are less steep and the 


number of points where the mode group velocity matches 
vp is reduced from 7 [cfr. Figs. 2(a) of the main text] to 
1. However, a large segment of the hybrid mode branch 
has approximately group velocity equal to vp in the range 
0.1 nm“^ ^ Q ^ 0-3 nm“^ and 0.10 eV < fici; < 0.15 eV. 
This results in an enhanced decay rate which appears as 
a broad peak of the imaginary part of the self-energy at 
hjjjeY ^ 1.5. 

The dependence of the decay rate (proportional to the 
imaginary part of the quasiparticle self-energy) on the 
Fermi energy is studied in detail in Fig. 4. We see that 
the decay rate depends smoothly on the Fermi energy 
and that there is no abrupt change in the imaginary 
part of the quasiparticle self-energy when the Fermi en¬ 
ergy crosses the extremes of the upper reststrahlen band. 
Moreover, for large Fermi energies, the main peak in the 
imaginary part of the self-energy is due to a plasmon- 
polariton mode which exists also in the absence of the 
hBN. This mode corresponds to the higher-energy branch 
in Fig. 2(a) of the main text and in Fig. 2 and generates 
the spectral feature labeled by (2) in Fig. 1 of the main 
text. To see this, in Fig. 4 we plot (horizontal dashed 
line) the expected position of the plasmon-polariton peak 
in the absence of the hBN slab. This is easily found 
from the well-known dispersion of the graphene plas¬ 
mon, which has group velocity equal to vp at the wave 
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FIG. 5. (Color online) Frequency dependence of the quasi¬ 
particle spectral function A{k,u}) for two values of k (indi¬ 
cated in the plots), corresponding to k = 0.01 and 0.5 kp- 
The Fermi energy is = 400 meV. The red solid lines corre¬ 
spond to d = 60 nm, while the black dashed line correspond 
to d = 0. Panel (a) The vertical green lines correspond to the 
expected positions of the plasmon-phonon polariton peaks for 
k 0. The two spectral features labeled (1) and (2) corre¬ 
spond to those in Fig. 1 of the main text. 


vector g* = j + Cb). Then, following the argu¬ 
ment explained in the main text, we find that the peak 
in the imaginary part of the quasiparticle self-energy is 
expected at hujje-F = —e^/{hvpe) + hvpq\\ — 1 —1.5, 

independent of other parameters. Satellite modes are 
clearly visible and drift to lower values of huu with de¬ 
creasing Fermi energy. The main peak broadens as the 
upper reststrahlen band is approached, and eventually a 
satellite peak becomes more prominent and closer to the 
frequency hwjsY ^ —1.5 of the main plasmon-polariton 
mode. 

Fig. 5 shows the frequency dependence of the quasi¬ 
particle spectral function at finite wave vector. Each 
panel shows the spectral function with (red solid line) 
and without (black dashed line) the hBN substrate. In 
Fig. 5(a), the quasiparticle peaks corresponding to con¬ 
duction and valence band, labeled by (1), are very close 
in energy and barely resolved. The spectral feature at 
hujjsp —1-7, labeled by (2), is due to the plasmon- 
polariton mode discussed above. The satellite peaks be- 
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FIG. 6. (Golor online) Difference between the spectral func¬ 
tion A(k,(ju) for d = 60 nm and the spectral function calcu¬ 
lated in the absence of the hBN slab, i.e. d = 0. The Fermi 
energy is ep = 400 meV. The dashed line indicates the Dirac 
cone. The redistribution of spectral weight due to the pres¬ 
ence of the hBN slab is sizable. 


tween (1) and (2) appear close to the frequencies where 
the peaks of the imaginary part of the quasiparticle self¬ 
energy are present. However, the expected position of the 
plasmon-polariton peak does not match its actual posi¬ 
tion (2). This effect is well-known and is a consequence of 
the attractive interaction between the emitted plasmon 
and the hole left behind by the photo-excited electron. [7] 
For the satellite peaks this effect is smaller as the fre¬ 
quency uj is increased. This can be understood from the 
plot of the imaginary part of the self-energy [Fig. 2(b) 
in the main text]: the plasmon-polariton peak has the 
largest decay rate, which obviously changes the position 
where the spectral function is maximal. The decay rate 
of the satellite peaks decreases as the frequency uj is in¬ 
creased. Inspecting Fig. 1 of the main text, we also note 
that k 0 OT k < kp are the wave vector ranges where 
the satellite peaks could be more clearly seen, because at 
intermediate values of k they are obfuscated as they cross 
the broad spectral features (1) and (2). Fig. 5(b) shows 
that at intermediate values of k the peak corresponding 
to the conduction band is clearly seen, while the valence 
band peak is strongly damped and disappears in the con¬ 
tinuum of excitations for huj ^ ep. 

Finally, Fig. 6 gives a more complete representation 
of the data shown in Fig. 5, by emphasizing the differ¬ 
ence between the quasiparticle spectral function with and 
without the hBN substrate. The effect of the substrate is 
sizable, with variations in magnitude of the spectral func¬ 
tion of order eV/ h. The spectral weight of both peaks 
(1) and (2) moves slightly to lower energies (i.e. the blue 
regions “move” to the red regions when d has a finite 
thickness). Moreover, the fine satellite bands are very 
clearly seen as an alternation of red and blue regions in 
this plot. 
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